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ABSTRACT 

We propose two efficient numerical methods of evaluating the luminosity distance in 
the spatially flat ACDM universe. The first method is based on the Carlson symmetric 
form of elliptic integrals, which is highly accurate and can replace numerical quadra- 
tures. The second method, using a modified version of Hermite interpolation, is less 
accurate but involves only basic numerical operations and can be easily implemented. 
We compare our methods with other numerical approximation schemes and explore 
their respective features and limitations. Possible extensions of these methods to other 
cosmological models are also discussed. 
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1 INTRODUCTION 

The computation of cosmological distances naturally arises 
in the study of cosmology, for example the luminosity dis- 
tance cLl in the analysis of type la supernova (SNIa) data 
and the angular diameter distance <1a in the study of gravi- 
tational lensing. These distances depend on the underlying 
cosmological model and their parameters. Therefore they 
are useful as cosmological tests. As a result, accurate and 
efficient numerical algorithms of evaluating these distances 
become a necessity for the practitioners of cosmological re- 
search. 

The analytical form of the cosmological distances can 
be derived from the solution of the Friedmann equation, 
an ordinary differential equation involving the scale factor 
a as a function of cosmic time t. Therefore, formulae for 
the distances usually involve an integral over the expansion 
history expressed in terms of the redshift z and cosmological 
parameters. In general, the integrations can be evaluated 
numerically by quadrature algorithms. However, numerical 
quadratures tend to be computationally heavy when high 
accuracy is desired. 

In the presence of this performance issue, it is advan- 
tageous to develop algorithms that are restricted to specific 
cosmological models but are otherwise more efficient than 
general-purpose quadratures. For the spatially flat ACDM 
model, efficient alg orithms for the luminosity distnace have 
been proposed by Pen ( 199 ^, he nceforward Pen99) and 
IWickramasinghe fc Ukwattal 1 20101 . henceforward WU10). 

In this paper we propose two different numerical meth- 
ods for the luminosity distance also in the context of the 
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spatially flat ACDM universe. The methods are presented 
in Sections 2 and 3 respectively. In Section 4, the perfor- 
mances of these methods are discussed. Finally in Section 
5 we discuss some possible extensions to the methods pre- 
sented in this paper. Throughout the paper we will focus 
on the luminosity distance only, but the results can be triv- 
ially extended to compute the angular diameter distance 
d A = d L /(l + z) 2 . 



2 METHOD I: EVALUATION USING 
CARLSON SYMMETRIC FORMS 

The luminosity distance di, in the spatially flat ACDM uni- 
verse is given by 



c(l + z) 
Ho 



dt 



V fi m(l+t) 3 +«A' 



(1) 



where fi m and Ha are the energy densities corresponding 
to the matter and cosmological con stant r espectively: fi m + 
JIa = 1. Following the notation in |Pen99l we introduce the 
parameter s — — S7 m )/r2 m and the change-of-variable 
u = 1/t, and re- write equation (JXJ) as 

d L l + z 



C/H y'sfLn 

where 
T(x) = 



T{s) - T 



l + z 



du 



Vm 4 + u 



(2) 



(3) 



The integral in equations JT} and ((3| are special 
cases of elliptic integrals. All elliptic integrals can be re- 
duced to several basic forms, the best known of which 
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are probably the thre e kind s of Legendre elliptic integrals 
ijWhittaker fc Watsonl Il969l . Chapt er 22), with reduction 
theor ems and examples presented in (|Abramowitz fc Stegunl 
1 1971 Chapter 17). In our case it is clearer to express this in- 
tegral by one of the Carlson symmetric forms Rf(xi, x 2 , x 3 ), 
which is defined as 



Rf{x 1 ,x 2 ,x 3 ) 



dt 



(4) 



y/(t + Xl){t + X 2 ){t + X 3 ) 

Using the reduction theorem^, it is straightforward to verify 



that 



m + 



(5) 



where 



m(x) 



-x + 1 2 

X X 



It has been known that the Carlson forms c an be 
computed numerically with high accuracy. ICarlsonl (|l979l ) 
showed that the computation of Rf can be accomplished 
iteratively with the error decreasing by a factor of 4 6 after 
each iteration, therefore achieving fast convergence. Further 
analysis of the al gorithms for R f and other elliptic integrals 
can be found in ( Car lsonll 1994 ). and computer implementa - 
tion details have been discussed in (|Carlson fc Notislll98ll ) 
and (|Press et al.ll2007l . Chapter 6). 



3 METHOD II: APPROXIMATION BY A 
MODIFIED HERMITE INTERPOLATION 

The method presented in Section 2 uses an iterative ap- 
proach to the computation of Rf- However, there are situ- 
ations where a closed, approx imate f ormula for the integral 
in equation ((3| is desired. In |Pen99l an app roxima tion was 
obtained using polynomial fit for T(x). In IWUld another 
method with higher accuracy was proposed. In this section 
we show how a modified version of Hermite interpol ation ca n 
lead to a class of approximations similar to that in |Pen99l . 

We intend to approximate equation Q using only ele- 
mentary operations, such as polynomial evaluation and nth 
root where n is a small integer. We note that the behavior 
of T(x) has several deficiencies. First, the derivative of T(x) 
becomes singular as x — > + . Second, the domain of T(x) 
extends to infinity. Either one is detrimental to the approxi- 
mation using polynomials. However, they can be removed by 
certain change-of-variables. For example, we can introduce 
a new function 



£(z) = T 2 (±-l) 



(6) 



that has smooth derivatives within the interval < x < 1 
and can be extended to the cases of x — > + and x — > 1~ . 
The limiting behaviors of £(ar) are shown below: 



e'(o+) 



-2A 



*(!") = 



0, 



where A = T(+oo) = 2.80436 • ■ ■ is a numerical constant- 
Using the end-point conditions in equation ((7j) one can 
construct a 3rd-order polynomial, which is a linear combi- 
nation of the four Hermite basis splines in [0, 1] , as a crude 
approximation with ~20% relative error. This linear com- 
bination is unique, allowing no further improvements. How- 
ever, we note that for realistic values of ilm it is not necessary 
to approximate in the entire interval [0, 1], because the 
subinterval [0, j4rr) corresponds to the scenario of 2 < 0, 
i.e. "the future". Therefore, we can introduce a free param- 
eter x* as the alternative lower end-point, and only perform 
the approximation in the subinterval [a;*, 1], if a constraint 
is put on Qm (or equivalently, s). 

To accommodate further refinements, a correction term 
w(x) can be added to the Hermite approximation. We re- 
quire the value and first derivative of w(x) to vanish at either 
end-point, so that it can be added to the Hermite approxi- 
mation without altering the coefficients on the basis splines. 
One choice of w(x) is made possible by a family of functions 



x 2 {l- 



(ax + b + 2a) 



(8) 



where a and b are adjustable parameters accounting for the 
deviation of the Hermite approximation from the true func- 
tion. Other choices are possible, but we will begin with the 
simple case of equation ©. 

By construction, the approximation described above 
has the property that the approximating function coincides 
with the true function at the end-points, a;* and 1, up to the 
first derivative. But we note that the goal is to approximate 
equation (J2J) rather than equation (J3J). This suggests that the 
implicit requirement of the coincidence of function values at 
end-points could be unnecessarily strong. Alternatively, we 
may refrain from requiring the approximating function val- 
ues to match the true ones. Instead, we only require the 
matching of first derivatives at x — x,, and leave the end- 
point value at another free parameter. To summarize, we 
now have four free parameters that can be tuned: x„, a and 
b, and the function value vo at x = x,. The approximation 
to equation © can be expressed as 

i(y) = v H^(y) + l[d H^(y)-AH[ 1 \y)] +w(y), (9) 
where 



I 



d = f'(a:«), 



and flf 5 are the Hermite basis splines, 

H { °\y) = 2y 3 -3y 2 + l, 
H?\y) = -2y 3 + 3y 2 , 



y A - 2y 2 + y, 



H[ r> (y) = y 3 -y 2 . 

Following the approach in |Pen99l . we choose the objec- 
tive function as the maximum relative error in di using the 
approximation ([Eq. [9]), with the restriction 0.2 ^ f2 m ^ 1. 
Minimizing the objective function over the parameters, we 
obtain the best-fit x* = 0.40176, a = 1.62053, b = -6.34985, 



(7) 



2 We note in passing that the constant X in |Pen99l . equation (5) 

1 See jOlver et alj l20ld . Chapter 19), available online at is identical to 1/A. A typo was made therein, which should have 

|http : //dlmf . nist . gov/19 ■ 29 1 been X = [J °°du /Vu 4 + u] ~ 1 . 
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and vo = 4.64111. Substituting the numerical values into 
equation ((9|, we therefore construct an approximation poly- 
nomial 

i(y) = 1.62053?/ 5 - 6.34985y 4 + 8.41443y 3 

-2.01328y 2 - 6.31293y + 4.64111. (10) 

Equation (|10|l is the main result of this section. With the pa- 
rameters determined, the approximation to can be com- 
puted using this formula with equations ([2]) and ([6]). 



4 PERFORMANCE OF THE METHODS 

In this section we proceed to assess the performance of the 
methods in Sections 2 and 3. The assessment is mainly done 
in terms of the accuracy and efficiency. 

4.1 Accuracy 

The first method can be used to yield highly accurate nu- 
merical approximation of dt for vast ranges of z and the 
parameter s if we adopt the algorithm for Rf by ICarlsonl 
l|l979l , ll994h . Unlike the methods based on the evaluation of 
a closed approximation formula, the desired cutoff error can 
be prescribed to determine when the iterative computation 
of Rf terminates. In practice, we found that the prescription 
of relative error ~10 -16 can be achieved without suffering 
significant loss in the computation speed. 

For the second method, we plot the distribution of the 
relative error of d_L in Figure [T] As can be seen from the fig- 
ure, the second method remains an approximation at best. 
Under our choice of fitting parameters and range of f2 m , the 
relative error in d_t is ~0.5%. The major source of this error 
is contributed by z < 0.1. For .1 < z < 10 our method 
is comparable with that of |Pen99l and ours slightly outper- 
forms it when z is larger. 

4.2 Efficiency 

Theoretically, the best-, worst-, and average-case temporal 
efficiencies for each method can be calculated or estimated 
by tracking every operation taken during the course of the 
computing. However, such a thorough analysis is beyond 
the scope of this paper. Instead, we empirically compare 
the running time of t he com pute r progr ams using the two 
methods with those of |Pen99l and lWUld under a controlled 
environment. 

In Figure [2] we display th e bench ma rk resu lts of our 
methods compared with that of |Pen99l and lWU10l . To simu- 
late a "real-world" application of these methods, we creates 
a sampl e of SNIa redshifts using the Supernova/Acceleration 
Probe (|SNAP Collaboration! [20041 ) fiducial redshift distri- 
bution containing 19 98 redshift points dis tributed within 
0.1 < z «S 1.7 fsee IShafieloo et all [20061 . Table 1). Our 
sample satisfies the same distribution to the SNAP fidu- 
cial, but is 16 times as dense, i.e. with 31968 points in to- 
tal. We have made c ustom implementations of the methods 
from |Pen99l . IWUlOl and our Method II in the C program- 
ming language, and uses the GNU Scientific Library (GSlB) 

3 http : //www . gnu . org/sof twar e/gsl/ 
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Figure 1. Contour plot for the distribution of the relative error 
in <1l using the approximation method in Section 3. Positive and 
negative values of the error are plotted in solid and dashed lines 
respectively. The "peaks" and "pits" in the left region of this 
figure (z < 0.1) dominate the global error. 

imple ment at ion of the Rf algorithm in (jCarlson fc Notisl 
Il98lf ) for Method I. In our benchmark routine, the com- 
puting of d_L values from our redshift sample is performed 
for n m = 0.3, 0.5, 0.7, and0.9 respectively, with each pass 
through the z sample repeated for 25 times (that sums up 
to a total of 3.2 x 10 6 evaluations of dz,). The benchmark 
itself is repeated for 2400 times. 

To interpret Figure [5] we make two remarks. First, the 
execution time results were collected from the output of the 
gprof profileiQ and does not reflect the absolute time spent. 
It is only meaningful as a relative measure useful for com- 
paring the speed of the codes with each other. Second, the 
results are dependent on our particular implementations as 
well as the computing environment. This is e vident if our 
Figure [2] is compared with Figure 4 in IWUlOl that s hows a 
reve rsed re sult for the speeds of the two methods in |Pen99l 
and lWUld 



5 DISCUSSION 

As Figure[5]suggcsts, both methods proposed in this paper is 
slower than the |Pen99l method. However, Method I is a very 
reasonable trade-off between an enormous gain in accuracy 
and small loss of efficiency. With Method I one does not need 
to resort to the numerical quadrature for the same level of 
accuracy. 

Method I can be extended to cover the ACDM model 



4 http : //www. gnu. org/sof tware/binutils/ 
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Figure 2. Histograms of the benchmark results for the four meth- 
ods. Dotted, thick solid, shaded, and thin solid histograms rep- 
resent the execution timi ngs of c odes implementing |Pen99l , our 
Method I, Method II, and lWUld respectively. Each histogram is 
normalized so that the total probability (the area enclosed under 
the boundaries) is unity. 



with a curvature term S7k, because in that case the equiva- 
lent of equation |TJ assumes the form 



d L 



1 + z 



c/H Q 
where 



dt 



(11) 



E(z-,n m ,n k ) = ^n m (i + z y + n k (i + z y + n A 

is the expansion rate (fL\ = 1 — — f^k), and 

fsin(x), fik < 0; 
sinn(:r) = I x. fik = 0; 

I sinh(x), f2k > 0. 

The integral in equation (| 1 1 1) is also an elliptic integral and 
can be reduced to Rf accordingly. This is potentially useful 
for the analysis of future SNIa data, because it has been sug- 
gested that the spacetime curvature should not be ignored 
in the probe of dark energy using lumin osity distance data 
jClarkson et al.ll2007l ; I6ztas et al.ll200Sl ). 

In contrast, Method II may not be as promising, be- 
cause i n its cu rrent form the accuracy does not outperform 
that of |Pen99l However, the idea behind the method may 
be useful when extending to alternative cosmological mod- 
els (for example, dynamical dark energy) which may not 
be reduced to the applicable scenarios of Method I. In the 
description of this method we have left some arbitrariness 
unjustified, notably the particular choice of the singularity- 
removing transformation (Eq. [6]), the parameterization of 
the correction term (Eq. [8]), and the very choice of Hermite 
basis splines. Alternative choices of them may be adopted to 



generate better approximations, for instance, the use of low- 
order Hermite-Birkhoff interpolation^ to selectively choose 
the point x £ [0, 1] near which the derivative information 
of the true function is to be best preserved. Moreover, our 
M ethod I I uses only elementary numerical operations, while 
m Iwuid the numerical logarithm is extensively used. 
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